/*
 * graph.c
 *
 * this module provides functions to manipulate the adjacency matrix of a graph.
 *
 *  Created on: Aug 16, 2011
 *      Author: correa
 */

#include <string.h>
#include "graph.h"

#define declareListForAssignment(r,k) 	int _i_; \
							int _next_[r+1]; \
							for (_i_ = 0; _i_ <= r; _i_++) _next_[_i_] = r; \
							int _covhead_[k+1]; \
							for (_i_ = 0; _i_ <= k; _i_++) _covhead_[_i_] = r

#define callAssignColor(r, R, k, cover, i, coversize)	assignColorByMatrix(r, R, k, cover, i, &coversize, _next_, _covhead_)
#define callAssignClique(r, R, k, cover, i, coversize)	assignCliqueByMatrix(r, R, k, cover, i, &coversize, _next_, _covhead_)

// sort specified vertices of a specified graph in a non-increasing order of degrees, in the
// following sense:
// in a subgraph of the specified graph induced by a set of vertices V = {V_1,V_2, ... ,V_i}, for 1 ≤ i ≤ r,
// it holds that V_i always has the maximum rank and the minimum degree in V.
// exdeg is used to break ties
// 2 vectors of r integers are defined in the C function stack
// vertices R[j] having rank[j]=min { rank[i] | i=0,1,...,r-1 } induce a regular subgraph
// mcr order is obtained by sorting R in a non-increasing order of rank
int mcrSort(const adjMatrix * g, const int r, const int * d, const int ** ng, int * rank) {
	setAdjMatrix(g);

	// degrees used locally
	int locald[r];
	memcpy(locald, d, r*sizeof(int));
	// used to break ties
	int exdeg[r];
	// a negative number indicates that this parameter is not known
	memset(exdeg, -1, r*sizeof(int));

	// determine min deg
	int Rmin[r];
	int nrmin = 1;
	int i = r - 1;

	int imin = 0;
	Rmin[0] = imin;
	int degrmin = locald[imin];
	int j;
	for (j = imin + 1; j < r; j++)
		if (locald[j] < degrmin) {
			nrmin = 1;
			Rmin[0] = j;
			degrmin = locald[j];
		}
		else if (locald[j] == degrmin)
			Rmin[nrmin++] = j;

	int rmin;
	memset(rank, 0, r*sizeof(int));
	while (i >= nrmin) {
		rmin = 0;
		if (nrmin > 1) {
			for (j = 1; j < nrmin; j++) {
				if (exdeg[Rmin[rmin]] < 0) {
					exdeg[Rmin[rmin]]=0;
					int l;
					for (l = 0; l < locald[Rmin[rmin]]; l++)
						if (!rank[ng[Rmin[rmin]][l]])
							exdeg[Rmin[rmin]] += locald[ng[Rmin[rmin]][l]];
				}
				if (exdeg[Rmin[j]] < 0) {
					exdeg[Rmin[j]]=0;
					int l;
					for (l = 0; l < locald[Rmin[j]]; l++)
						if (!rank[ng[Rmin[j]][l]])
							exdeg[Rmin[j]] += locald[ng[Rmin[j]][l]];
				}
				if (exdeg[Rmin[j]] < exdeg[Rmin[rmin]])
					rmin = j;
			}
		}

		rank[Rmin[rmin]] = --i;

		while (rank[imin])
			imin++;
		for (j = 0; j < locald[Rmin[rmin]]; j++)
			if (!rank[ng[Rmin[rmin]][j]]) {
				exdeg[ng[Rmin[rmin]][j]] -= locald[Rmin[rmin]];
				locald[ng[Rmin[rmin]][j]]--;
			}

		Rmin[0] = imin;
		nrmin = 1;
		degrmin = locald[imin];
		for (j = imin + 1; j < r; j++)
			if (!rank[j]){
				if (locald[j] < degrmin) {
					nrmin = 1;
					Rmin[0] = j;
					degrmin = locald[j];
				}
				else{
					if (locald[j] == degrmin)
						Rmin[nrmin++] = j;
				}
			}
	}

	--i;
	for (j = imin; j < r; j++)
		if (!rank[j])
			rank[j] = i;

	return 0;
}

// sort specified vertices of a specified graph in a non-increasing order of degrees, in the
// following sense:
// in a subgraph of the specified graph induced by a set of vertices V = {V_1,V_2, ... ,V_i}, for 1 ≤ i ≤ r,
// it holds that V_i always has the maximum rank and the minimum degree in V.
// exdeg is used to break ties; maximum rank is returned
// 3 vectors of r integers are defined in the C function stack
// if i is returned, then vertices R[j] having rank[j]=i induce a regular subgraph
int mcrDegreeSort(const adjMatrix * g, const int r, const int * R, int * rank) {
	setAdjMatrix(g);

	// degrees used locally
	int i, j;
	int deg[r];
	memset(deg, 0, r*sizeof(int));
	for (i = 0; i < r; i++) {
		for (j = i + 1; j < r; j++)
			if (!hasEdge(R[i], R[j])) {
				deg[i]++;
				deg[j]++;
			}
	}
	// used to break ties
	int exdeg[r];
	// a negative number indicates that this parameter is not known
	memset(exdeg, -1, r*sizeof(int));

	// determine min deg
	int Rmin[r];
	int nrmin = 1;
	i = r - 1;

	int imin = 0;
	Rmin[0] = imin;
	int degrmin = deg[imin];
	for (j = imin + 1; j < r; j++)
		if (deg[j] < degrmin) {
			nrmin = 1;
			Rmin[0] = j;
			degrmin = deg[j];
		}
		else if (deg[j] == degrmin)
			Rmin[nrmin++] = j;

	int rmin;
	memset(rank, 0, r*sizeof(int));
	while (i >= nrmin) {
		rmin = 0;
		if (nrmin > 1) {
			for (j = 1; j < nrmin; j++) {
				if (exdeg[Rmin[rmin]] < 0) {
					exdeg[Rmin[rmin]]=0;
					int l;
					for (l = imin; l < r; l++)
						if (!rank[l] && l != Rmin[rmin] && hasEdge(R[Rmin[rmin]], R[l]))
							exdeg[Rmin[rmin]] += deg[l];
				}
				if (exdeg[Rmin[j]] < 0) {
					exdeg[Rmin[j]]=0;
					int l;
					for (l = imin; l < r; l++)
						if (!rank[l] && l != Rmin[j] && hasEdge(R[Rmin[j]], R[l]))
							exdeg[Rmin[j]] += deg[l];
				}
				if (exdeg[Rmin[j]] < exdeg[Rmin[rmin]])
					rmin = j;
			}
		}

		rank[Rmin[rmin]] = --i;

		while (rank[imin])
			imin++;
		for (j = imin; j < r; j++)
			if (!rank[j] && hasEdge(R[Rmin[rmin]], R[j])) {
				exdeg[j] -= deg[Rmin[rmin]];
				deg[j]--;
			}

		Rmin[0] = imin;
		nrmin = 1;
		degrmin = deg[imin];
		for (j = imin + 1; j < r; j++)
			if (!rank[j]){
				if (deg[j] < degrmin) {
					nrmin = 1;
					Rmin[0] = j;
					degrmin = deg[j];
				}
				else if (deg[j] == degrmin){
					Rmin[nrmin++] = j;
				}
			}
	}

	--i;
	for (j = imin; j < r; j++)
		if (!rank[j])
			rank[j] = i;

	return i;
}

// sort specified vertices of a specified graph in a non-increasing order of anti-degrees, in the
// following sense:
// in a subgraph of the specified graph induced by a set of vertices V = {V_1,V_2, ... ,V_i}, for 1 ≤ i ≤ r,
// it holds that V_i always has the maximum rank and the minimum anti-degree in V.
// exantideg is used to break ties; maximum rank is returned
// 3 vectors of r integers are defined in the C function stack
// if i is returned, then vertices R[j] having rank[j]=i induce a regular subgraph
int mcrAntiDegreeSort(const adjMatrix * g, const int r, const int * R, int * rank) {
	setAdjMatrix(g);

	// anti-degrees used locally
	int i, j;
	int antideg[r];
	memset(antideg, 0, r*sizeof(int));
	for (i = 0; i < r; i++) {
		for (j = i + 1; j < r; j++)
			if (!hasEdge(R[i], R[j])) {
				antideg[i]++;
				antideg[j]++;
			}
	}
	// used to break ties
	int exdeg[r];
	// a negative number indicates that this parameter is not known
	memset(exdeg, -1, r*sizeof(int));

	// determine min antideg
	int Rmin[r];
	int nrmin = 1;
	i = r - 1;

	int imin = 0;
	Rmin[0] = imin;
	int antidegrmin = antideg[imin];
	for (j = imin + 1; j < r; j++)
		if (antideg[j] < antidegrmin) {
			nrmin = 1;
			Rmin[0] = j;
			antidegrmin = antideg[j];
		}
		else if (antideg[j] == antidegrmin)
			Rmin[nrmin++] = j;

	int rmin;
	memset(rank, 0, r*sizeof(int));
	while (i >= nrmin) {
		rmin = 0;
		if (nrmin > 1) {
			for (j = 1; j < nrmin; j++) {
				if (exdeg[Rmin[rmin]] < 0) {
					exdeg[Rmin[rmin]]=0;
					int l;
					for (l = imin; l < r; l++)
						if (!rank[l] && l != Rmin[rmin] && !hasEdge(R[Rmin[rmin]], R[l]))
							exdeg[Rmin[rmin]] += antideg[l];
				}
				if (exdeg[Rmin[j]] < 0) {
					exdeg[Rmin[j]]=0;
					int l;
					for (l = imin; l < r; l++)
						if (!rank[l] && l != Rmin[j] && !hasEdge(R[Rmin[j]], R[l]))
							exdeg[Rmin[j]] += antideg[l];
				}
				if (exdeg[Rmin[j]] < exdeg[Rmin[rmin]])
					rmin = j;
			}
		}

		rank[Rmin[rmin]] = --i;

		while (rank[imin])
			imin++;
		for (j = imin; j < r; j++)
			if (!rank[j] && hasEdge(R[Rmin[rmin]], R[j])) {
				exdeg[j] -= antideg[Rmin[rmin]];
				antideg[j]--;
			}

		Rmin[0] = imin;
		nrmin = 1;
		antidegrmin = antideg[imin];
		for (j = imin + 1; j < r; j++)
			if (!rank[j]){
				if (antideg[j] < antidegrmin) {
					nrmin = 1;
					Rmin[0] = j;
					antidegrmin = antideg[j];
				}
				else if (antideg[j] == antidegrmin){
					Rmin[nrmin++] = j;
				}
			}
	}

	--i;
	for (j = imin; j < r; j++)
		if (!rank[j])
			rank[j] = i;

	return i;
}

static inline int assignByList(const int * d, const int ** ng, const int k, int * cover, const int i, int * coversize, int * clique) {
	int j;
	// current clique starts with size 1
	cover[i] = 0;

	// find clique to include vertex as the smallest
	// not appearing in the neighborhood of imaxantideg
	memset(clique, 0, ((*coversize)+1)*sizeof(int));
	for (j = 0; j < d[i]; j++)
		clique[cover[ng[i][j]]] = 1;

	int c = 1;
	while (c <= *coversize && clique[c])
		c++;
	if (c <= k) {
		if (c > *coversize)
			(*coversize)++;
		cover[i] = c;
	}

	return c;
}

int greedySort(const adjMatrix * g, const int r, const int * d, const int ** ng, const int k, int * cover) {
	setAdjMatrix(g);

	int coversize = 0;
	int i;
	 // array, indexed by the cliques, used to mark all cliques in the neighborhood of a vertex
	int clique[k+1];

	memset(cover, 0, r*sizeof(int));
	for (i = 0; i < r; i++) {
		assignByList(d, ng, k, cover, i, &coversize, clique);
	}

	return coversize;
}

static inline int assignColorByMatrix(const int r, const int * R, const int k, int * cover, const int i, int * coversize, int * next, int * covhead) {
	int c = 0;
	int cc = covhead[c];
	int antcc = r;
	int vai = R[i];
	while (!cover[i] && c < *coversize){
		// search for cover to add
		for (; cc != covhead[c+1] && !hasEdge(vai, R[cc]); antcc = cc, cc = next[cc]);
		if (cc == covhead[c+1])
			cover[i] = c + 1;
		else {
			c++;
			for (; cc != covhead[c]; antcc = cc, cc = next[cc]);
		}
	}
	if (c == *coversize && *coversize < k) {
		(*coversize)++;
		covhead[c] = i;
		cover[i] = c + 1;
	}
	// include i in the tail of cover list
	next[antcc] = i;
	next[i] = cc;

	return c + 1;
}

static inline int assignCliqueByMatrix(const int r, const int * R, const int k, int * cover, const int i, int * coversize, int * next, int * covhead) {
	int c = 0;
	int cc = covhead[c];
	int antcc = r;
	int vai = R[i];
	while (!cover[i] && c < *coversize){
		// search for cover to add
		for (; cc != covhead[c+1] && hasEdge(vai, R[cc]); antcc = cc, cc = next[cc]);
		if (cc == covhead[c+1])
			cover[i] = c + 1;
		else {
			c++;
			for (; cc != covhead[c]; antcc = cc, cc = next[cc]);
		}
	}
	if (c == *coversize && *coversize < k) {
		(*coversize)++;
		covhead[c] = i;
		cover[i] = c + 1;
	}
	// include i in the tail of cover list
	next[antcc] = i;
	next[i] = cc;

	return c + 1;
}

// k: maximum number of sets
// 1 vector of r integers are defined in the C function stack
// 1 vector of k integers are defined in the C function stack
int greedyColoringSort(const adjMatrix * g, const int r, const int * R, const int k, int * cover) {
	setAdjMatrix(g);

	int coversize = 0;
	declareListForAssignment(r, k);

	int i;
	for (i = 0; i < r; i++) {
		cover[i] = 0;
		callAssignColor(r, R, k, cover, i, coversize);
	}

	return coversize;
}

// sort R according to a clique cover of the subgraph induced by R
// clique cover is determined using greedy algorithm in the vertex order specified by the input value of R
// 1 vector of r integers are defined in the C function stack
// it returns the size of the clique cover and sets the cover array such that:
// for i in 0, ..., r-1, the entry cover[i] contains the size of a clique
// cover of the subgraph induced by R[0], ..., R[i]
int greedyCliqueCoverSort(const adjMatrix * g, const int r, const int * R, const int k, int * cover) {
	setAdjMatrix(g);

	int coversize = 0;
	declareListForAssignment(r, k);

	int i;
	for (i = 0; i < r; i++) {
		cover[i] = 0;
		callAssignClique(r, R, k, cover, i, coversize);
	}

	return coversize;
}

// sort R according to a clique cover or coloring of the subgraph induced by R,
// depending whether d and ng represents neighborhoods or antineighborhoods.
// index in R is used as last break tie condition
// Assuming they represent antineighborhoods, clique cover is determined using dsatur
// 4 vectors of r integers are defined in the C function stack
// it returns the size of the clique cover and sets the cover array such that:
// for i in 0, ..., r-1, the entry cover[i] contains the size of a clique
// cover of the subgraph induced by all vertices of rank less than or equal to cover[i]
// the number of cliques in the cover is returned. Cliques returned in the cover array are numbered from 1 to number of cliques in the
int dsaturSort(const adjMatrix * g, const int r, const int * d, const int ** ng, const int k, int * cover) {
	setAdjMatrix(g);

	int next[r+1];
	int prev[r+1];
	int satdeglist[r];
	int satdeg[r];
	// array, indexed by the cliques, used to mark all cliques in the neighborhood of imaxantideg
	int clique[k+1];
	int nclique = 0;

	int i, j;

	// initially, all cliques are empty, all sauration degrees are 0 and
	// its doubly linked list contains all vertices
	memset(cover, 0, r*sizeof(int));
	memset(satdeg, 0, r*sizeof(int));
	for (i = 0; i < r; i++) {
		next[i] = i+1;
		prev[i] = i-1;
		satdeglist[i] = r;
	}
	satdeglist[0] = 0;
	prev[0] = r;

	int ncovered;
	int maxsatdeg = 0;
	int imaxantideg;
	for (ncovered = 0; ncovered < r; ncovered++) {
		// find vertex with max antideg among those with max satdeg
		i = satdeglist[maxsatdeg];
		imaxantideg = i;
		while (next[i] < r) {
			i = next[i];
			if (d[i] > d[imaxantideg])
				imaxantideg = i;
		}

		// remove vertex form its doubly linked list
		if (prev[imaxantideg] == r) {
			satdeglist[maxsatdeg] = next[imaxantideg];
			if (ncovered < r - 1)
				while (satdeglist[maxsatdeg] == r)
					maxsatdeg--;
		}
		else
			next[prev[imaxantideg]] = next[imaxantideg];
		prev[next[imaxantideg]] = prev[imaxantideg];

		// find clique to include vertex as the smallest
		// not appearing in the neighborhood of imaxantideg
		int c = assignByList(d, ng, k, cover, imaxantideg, &nclique, clique);

		// update satdeglist
		if (c <= k)
			for (i = 0; i < d[imaxantideg]; i++) {
				int ii = ng[imaxantideg][i];
				if (!cover[ii]) {
					for (j = 0; j < d[ii] && (ng[ii][j] == imaxantideg || cover[ng[ii][j]] != c); j++);
					if (j == d[ii]) {
						// remove i from list
						if (prev[ii] == r)
							satdeglist[satdeg[ii]] = next[ii];
						else
							next[prev[ii]] = next[ii];
						prev[next[ii]] = prev[ii];

						// insert i in next list
						satdeg[ii]++;
						next[ii] = satdeglist[satdeg[ii]];
						prev[ii] = r;
						prev[satdeglist[satdeg[ii]]] = ii;
						satdeglist[satdeg[ii]] = ii;

						if (satdeg[ii] > maxsatdeg)
							maxsatdeg = satdeg[ii];
					}
				}
			}
	}

	return nclique;
}

// sort R according to a clique cover of the subgraph induced by R,
// using the symmetric adjacency matrix
// clique cover is determined using dsatur
// 5 vectors of r integers are defined in the C function stack
// it returns the size of the clique cover and sets the cover array such that:
// for i in 0, ..., r-1, the entry cover[i] contains the size of a clique
// cover of the subgraph induced by all vertices of rank less than or equal to cover[i]
// the number of cliques in the cover is returned. Cliques returned in the cover array are numbered from 1 to number of cliques in the
int dsaturCliqueCoverSort(const adjMatrix * g, const int r, const int * R, const int k, int * cover) {
	setAdjMatrix(g);

	declareListForAssignment(r, k);

	int antideg[r];
	int next[r+1];
	int prev[r+1];
	int satdeglist[r];
	int satdeg[r];
	int nclique = 0;

	int i, j;

	// initially, all cliques are empty, all saturation degrees are 0 and
	// its doubly linked list contains all vertices
	memset(cover, 0, r*sizeof(int));
	memset(satdeg, 0, r*sizeof(int));
	for (i = 0; i < r; i++) {
		next[i] = i+1;
		prev[i] = i-1;
		satdeglist[i] = r;
		for (j = i + 1; j < r; j++)
			if (!hasEdge(R[i], R[j])) {
				antideg[i]++;
				antideg[j]++;
			}
	}
	satdeglist[0] = 0;
	prev[0] = r;

	int ncovered;
	int maxsatdeg = 0;
	int imaxantideg;
	for (ncovered = 0; ncovered < r; ncovered++) {
		// find vertex with max antideg among those with max satdeg
		i = satdeglist[maxsatdeg];
		imaxantideg = i;
		while (next[i] < r) {
			i = next[i];
			if (antideg[i] > antideg[imaxantideg])
				imaxantideg = i;
		}

		// remove vertex form its doubly linked list
		if (prev[imaxantideg] == r) {
			satdeglist[maxsatdeg] = next[imaxantideg];
			if (ncovered < r - 1)
				while (satdeglist[maxsatdeg] == r)
					maxsatdeg--;
		}
		else
			next[prev[imaxantideg]] = next[imaxantideg];
		prev[next[imaxantideg]] = prev[imaxantideg];

		int c = callAssignClique(r, R, k, cover, imaxantideg, nclique);

		// update satdeglist
		if (c <= k)
			for (i = 0; i < r; i++)
				if (!hasEdge(R[imaxantideg], R[i]) && !cover[i]) {
					for (j = 0; j < r && !hasEdge(R[i], R[j]) && (j == imaxantideg || cover[j] != c); j++);
					if (j == r) {
						// remove i from list
						if (prev[i] == r)
							satdeglist[satdeg[i]] = next[i];
						else
							next[prev[i]] = next[i];
						prev[next[i]] = prev[i];

						// insert i in next list
						satdeg[i]++;
						next[i] = satdeglist[satdeg[i]];
						prev[i] = r;
						prev[satdeglist[satdeg[i]]] = i;
						satdeglist[satdeg[i]] = i;

						if (satdeg[i] > maxsatdeg)
							maxsatdeg = satdeg[i];
					}
				}
	}

	return nclique;
}

// sort R according to a coloring of the subgraph induced by R,
// using the symmetric adjacency matrix
// coloring is determined using dsatur
// 5 vectors of r integers are defined in the C function stack
// it returns the size of the coloring and sets the cover array such that:
// for i in 0, ..., r-1, the entry cover[i] contains the size of a coloring
// of the subgraph induced by all vertices of rank less than or equal to cover[i]
// the number of colors in the coloring is returned. Colors returned in the cover array are numbered from 1 to number of cliques in the
int dsaturColoringSort(const adjMatrix * g, const int r, const int * R, const int k, int * cover) {
	setAdjMatrix(g);

	declareListForAssignment(r, k);

	int deg[r];
	int next[r+1];
	int prev[r+1];
	int satdeglist[r];
	int satdeg[r];
	int ncolor = 0;

	int i, j;

	// initially, all cliques are empty, all sauration degrees are 0 and
	// its doubly linked list contains all vertices
	memset(cover, 0, r*sizeof(int));
	memset(satdeg, 0, r*sizeof(int));
	for (i = 0; i < r; i++) {
		next[i] = i+1;
		prev[i] = i-1;
		satdeglist[i] = r;
		for (j = i + 1; j < r; j++)
			if (hasEdge(R[i], R[j])) {
				deg[i]++;
				deg[j]++;
			}
	}
	satdeglist[0] = 0;
	prev[0] = r;

	int ncovered;
	int maxsatdeg = 0;
	int imaxdeg;
	for (ncovered = 0; ncovered < r; ncovered++) {
		// find vertex with max antideg among those with max satdeg
		i = satdeglist[maxsatdeg];
		imaxdeg = i;
		while (next[i] < r) {
			i = next[i];
			if (deg[i] > deg[imaxdeg])
				imaxdeg = i;
		}

		// remove vertex form its doubly linked list
		if (prev[imaxdeg] == r) {
			satdeglist[maxsatdeg] = next[imaxdeg];
			if (ncovered < r - 1)
				while (satdeglist[maxsatdeg] == r)
					maxsatdeg--;
		}
		else
			next[prev[imaxdeg]] = next[imaxdeg];
		prev[next[imaxdeg]] = prev[imaxdeg];

		int c = callAssignColor(r, R, k, cover, imaxdeg, ncolor);

		// update satdeglist
		if (c <= k)
			for (i = 0; i < r; i++)
				if (hasEdge(R[imaxdeg], R[i]) && !cover[i]) {
					for (j = 0; j < r && hasEdge(R[i], R[j]) && (j == imaxdeg || cover[j] != c); j++);
					if (j == r) {
						// remove i from list
						if (prev[i] == r)
							satdeglist[satdeg[i]] = next[i];
						else
							next[prev[i]] = next[i];
						prev[next[i]] = prev[i];

						// insert i in next list
						satdeg[i]++;
						next[i] = satdeglist[satdeg[i]];
						prev[i] = r;
						prev[satdeglist[satdeg[i]]] = i;
						satdeglist[satdeg[i]] = i;

						if (satdeg[i] > maxsatdeg)
							maxsatdeg = satdeg[i];
					}
				}
	}

	return ncolor;
}
